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1.  Technical  Summary  of  Research  Accomplished. 

We  will  first  briefly  describe  the  research  results  that  have  already  been  published  [1]  in  the  scientific 
literature,  and  then  we  will  provide  some  details  of  the  remaining  publications  that  resulted  from  the  research 
effort  and  have  been  submitted  for  publication  or  are  still  in  preparation  to  be  submitted. 

1.1.  Summary  of  Published  Research  Accomplishments. 

1.1.1.  Introduction.  The  propagation  of  short-duration  electromagnetic  pulses  in  dielectric  media  is 
of  interest  in  diverse  technological  applications  areas  such  as  geophysics  (ground-penetrating  radar  [2])  and 
bioelectromagnetics  (health  and  safety  analyses  [3]).  Most  importantly,  time-domain  reflectometry  (TDR, 
[4])  employs  short-duration  electromagnetic  pulses  in  a  laboratory  setting  to  obtain  the  frequency-dependent 
dielectric  properties  of  materials.  Time-domain  waveforms  of  propagated  signals  in  TDR  experiments  are 
routinely  measured  and  analyzed  [5]- [6].  In  the  TDR  setting,  the  experimentally  determined  complex¬ 
valued  dielectric  permittivity  e(cj)  is  subsequently  fitted  to  a  model  assumed  to  be  relevant  for  the  range  of 
frequencies  in  which  measurements  are  obtained.  The  two  models  commonly  used  to  describe  the  variation  of 
the  permittivity  with  frequency  up  to  about  10  GHz  are  the  Debye  [7]  and  the  Cole-Cole  [8]  models.  Examples 
of  the  use  of  the  Debye  and  Cole-Cole  models  to  fit  experimentally  determined  dielectric  permittivity  data  for 
human  tissue  include  [9]  and  [10]  respectively.  The  reasons  for  choosing  a  particular  model  to  fit  a  particular 
set  of  data  are  not  always  clear,  e.g.,  both  of  the  aforementioned  models  have  been  used  to  fit  permittivity 
data  obtained  for  human  tissue  [10]  and  concrete  [11].  A  modeled  small-  and  large-depth  response  of  an 
experimentally  examined  dielectric  subjected  to  any  physically  realizable  pulse  (hence  measurable  in  a  TDR 
experiment)  can  be  obtained  by  convolution  with  the  appropriate,  theoretically  obtained,  impulse  response 
function;  subsequent  comparisons  to  actual  response  waveforms  measured  in  a  TDR  setup  could  be  used  to 
determine  which  model  is  most  appropriate  to  represent  a  given  set  of  measured  permittivity  data. 

The  response  of  the  Debye  dielectric  medium  model  to  transient  electromagnetic  radiation  has  been 
investigated  in  [12]-[14] .  These  authors  found  that  the  wavefront  impulse  response  supports  discontinuities 
on  x  =  Coot,  where  c^  =  l/y/e(oo)fio  is  the  wavefront  speed,  and  that  it  decays  exponentially  with  depth 
hence  it  is  important  only  up  to  a  distance  of  O(coot)  past  the  air /dielectric  interface  (the  so-called  ”  time- 
domain  skin-depth”).  Also  [12]-[14],  it  was  found  that  past  that  thin  layer  the  impulse  response  satisfies  an 
advection-diffusion  equation,  that  it  travels  with  the  zero- frequency  phase  speed  along  the  sub-characteristic 
x  =  =  cot ,  where  es  =  6(0)  is  the  DC  permittivity,  and  that  it  diffuses  around  that  sub-characteristic. 

In  this  paper  we  examine  the  small-  and  large-depth  response  of  a  Cole-Cole  dielectric  half-space  sub¬ 
jected  to  a  prescribed  incident  pulse;  the  case  of  delta-function  incidence  is  employed  to  determine  and 
analyze  the  resulting  impulse  response.  Our  purpose  is  to  contrast  our  findings  to  the  corresponding  ones 
obtained  for  the  Debye  model  in  order  to  ascertain  whether  the  time-domain  waveforms  obtained  in  a  TDR 
experiment  could  serve  as  a  means  for  selecting  the  most  appropriate  frequency-domain  model  for  the  ex¬ 
perimentally  obtained  dielectric  data.  Our  approach  involves  both  asymptotic  and  numerical  methods.  We 
find  that  the  Cole-Cole  model’s  impulse  response  is  infinitely  smooth  at  the  wavefront  (small-depth),  and 
determine  its  shape.  It  follows  that  sawtooth  and  square-pulse  waveforms,  and  all  other  realistic  waveforms, 
become  smooth  after  travelling  a  brief  time  in  any  Cole-Cole  model.  This  is  in  contrast  to  the  case  of  the 
Debye  impulse  response  which  is  discontinuous  at  the  wavefront.  Also,  we  find  that  the  location  of  the  peak 
of  the  main  (large-depth)  response  in  the  Cole-Cole  dielectric  model  occurs  at  an  earlier  space-time  location 
than  that  found  for  the  Debye  dielectric  model  main  response. 

1.1.2.  Problem  Formulation.  The  dielectric  dispersion  proposed  in  [8]  defines  what  is  known  in  the 
literature  as  the  frequency-domain  Cole-Cole  dielectric  medium  model  which,  in  the  context  of  the  Laplace 
transform  (s  =  iuj ),  takes  the  form 


t(s)  =  6oo  + 


1  + (sr)a  ’ 


(1.1) 
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where  0  <  a  <  1  is  a  data- fit  parameter,  and  e^,  es,  r  are  respectively  the  infinite-frequency  permittivity,  the 
zero-frequency  permittivity,  and  the  central  relaxation  time  [15]  around  which  a  distribution  of  relaxation 
times  exists  and  whose  width  is  controlled  by  a.  The  four  parameters  in  (1.1)  are  fitted  to  dielectric 
permittivity  data  measured  over  a  frequency  band,  e.g.,  [10]  fits  frequency-domain  dielectric  permittivity 
data  obtained  for  various  biological  tissue  types  to  a  linear  combination  of  4  Cole-Cole  models.  It  is  known 
[16]  that  (1.1)  constitutes  a  causal  dielectric  model. 

For  simplicity,  we  will  consider  one-dimensional  electromagnetic  pulse  propagation  in  a  homogeneous, 
non-magnetic  (/i  =  /zo,  Mo  is  the  permeability  of  vacuum)  half-space,  x  >  0,  whose  frequency-dependent 
dielectric  permittivity  is  modeled  by  (1.1).  The  electric  field  E(x,  t)  is  prescribed  at  x  =  0.  The  time- 
dependent  Maxwell  system  reduces  to 


(1.2) 


d2E  d2p  i  d2E 
dt 2  dt2  Mo  dx2 


x  >  0,  t  >  0, 


where  P{x)  t)  is  the  induced  polarization  field  which  describes  the  response  of  the  dielectric  medium  to  the 
propagating  wave.  Using  (1.1),  the  induced  polarization  is  defined  as  [17] 


(1.3) 


Xa{t  —  t  )E{x ,  t  )dt  , 


t  >  0, 


where 

(1.4) 


Xn(t)  =C  1{ 


-<oo  .  =  j_  r+o° 
1  +  (sr)“  '  2? n  7c_ioo 


£s  ^OO 

1  + (sr)a 


estds , 


t  >  0, 


is  the  Cole-Cole  time-domain  susceptibility  kernel.  When  a  =  1,  the  Debye  time-domain  susceptibility  is 
obtained.  In  (1.4),  C~l  denotes  the  inverse  Laplace  transform  performed  along  the  standard  Bromwich 
contour  (7m{C}  =  0,  C  >  0)  on  the  complex  s-plane  cut  along  the  Rc{s}  <  0  axis  with  the  branch  point 
at  the  origin  so  that  —  n  <  arg(s)  <  7r,  i.e.,  the  denominator  does  not  introduce  any  additional  singularities 
in  the  complex  s  plane.  The  initial  value  P(x,  0)  =  0  is  evident  from  (1.3).  The  polarization  thusly  defined 
also  satisfies  the  fractional-order  differential  equation 


(1.5)  Ta^-  +  P  =  (es-eQO)E,  t>  0;  P(x,0)=0, 

i.e.,  Xa(t)  is  the  Green’s  function  for  (1.5)  which  is  obtained  when  E  =  6(t).  The  nonlocal-in-time  operator 
0  <  a  <  1,  denotes  the  Caputo  fractional-order  time  derivative  [18]  which,  for  functions  satisfying 
P(x,  0)  =  0,  is  defined  as 


daP(x,t)  _  1  r*  P\x,u)du 

1  J  dta  ~  r(l  -  a)  Jo  ( t  - «)“  ’ 

where  P  represents  differentiation  with  respect  to  the  time  argument  of  P.  Finally,  we  assume  that  at  t  =  0 
all  fields  are  zero  in  the  half-space  x  >  0,  that  all  fields  vanish  as  x  — >  oo,  t  >  0,  and  that  the  electric  field 
is  prescribed  at  x  =  0,  i.e.,  we  consider  a  signaling  problem  with  boundary  data 

(1.7)  £(0,  *)  =  /(*);  *>o, 

where  f(t)  is  the  prescribed  electric  field. 

1.1.3.  Asymptotic  Analysis  for  Pulse  Propagation.  We  use  the  Laplace  transform  in  time  to 
solve  (1.2)  and  (1.5)  subject  to  (1.7),  the  initial  conditions,  and  the  behavior  at  infinity  described  in  Section 
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2.  After  eliminating  the  P  field  we  determine  that  for  the  signaling  problem  in  a  Cole- Cole  dielectric  half 
space  (x  >  0)  the  electric  field  is  given  by  the  following  Bromwich  integral 


(1.8) 


E{x 


j  K+too 

’  ^  “  2m  7(_<00 


F(s)eS^  c<~  V *a+0-y}(is,  t  > 


where  -7r  <  arg(s)  <  7r,  Re.yJ  >  0,  0  =  7  =  and  F(s)  =  £{/(£)}.  Closing  the  Bromwich 

contour  with  a  semi-circle  to  the  right  of  s  =  £  we  obtain  E(x)  t)  =  0,  t  <  as  the  integrand  exhibits  no 
singularities  in  that  region  (causality).  For  the  branch  chosen,  the  square  root  in  (1.8)  does  not  introduce 
additional  singularities  (recall  0  <  a  <  1)  since  the  argument  of  the  roots  of  the  numerator  and  denominator 
is  arg(s)  =  (n  +  2/c7t)/q;  k  =  0,  ±1, i.e.,  these  roots  are  outside  the  chosen  principal  branch. 

We  next  consider  (1.8)  for  t  — ►  ,  i.e.,  for  times  immediately  after  the  arrival  of  the  wavefront  described 

by  the  characteristic  ray  x  =  Coot  at  a  fixed  x  location,  and  for  t  ^  3>  1. 

The  Wavefront  Region  (t  — j ► 

A  large— s  expansion  of  the  bracketed  expression  in  (1.8) 


(1.9) 

results  in 

(1.10) 


Coc  ^oo 


Sl±L)  =  (1  _  JLj  _  JL 

S“+/V  ^  Coo’  Coo  2 


K+ioo 

E(x,t)  =  /  F(s)e~As  "“(1  +  0(s1-2Q))e3(‘-^)ds, 

2m 


where  A  =  _2l  >  0.  Due  to0<l-a<l,we  determine  that 

Cqo  ^ 

dn 

lim  - — E(x,t)  =  0,  n  >  0, 

+  ,  1+  dtn 


since 


lim  5nF(s)(l  -hO(51"jQ))e 


=  0,  n  >  0. 


Thus,  the  impulse  response  of  the  Cole-Cole  medium  is  infinitely  smooth  at  the  wavefront  even  in  the  case 
of  F(s)  =  1  ( f(t )  =  6(t)  signaling  data).  In  contrast,  setting  a  =  1  (Debye  model)  in  (1.10)  we  obtain 

(1.11) 


E(x,t)&e  Vf(t — — ), 


i.e.,  the  Debye  wavefront  impulse  response  inherits  the  continuity  properties  of  /(£)  at  t  =  0  [13]. 

For  the  special  case  a  =  1/2,  standard  Laplace  transform  tables  allow  us  to  invert  the  leading-order 
term  in  (1.10)  via 


(1.12) 


,-l  1  0(1-7)  X 


C  l{e 


}  = 


4a/7T  C( 

and  obtain  the  approximate  short-time  response 

1  0(1-7) 


t 3/2 


Coot  ,  t  >  0, 


(1.13) 


E(x,  t) 


4 \fF 


Jx/Coo  U  X/i 


Coo)3/2 


^2(i-7)2x2 

iec2  _x/Coo) 


(a) 


(b) 


Fig.  1.1.  Wavefront  response  for  8(t)  signaling  data  and  0.5  <  or  <  0.8;  T  =  t  —  -2— . 


X=0.01  X*0.1 


(a) 


(b) 


Fig.  1.2.  The  wavefront  impulse  response  (1.15)  at  two  depths  for  a  =  0.7. 


For  other  values  of  a  we  proceed  by  collapsing  the  Bromwich  contour  in  (1.10)  onto  the  branch  cut 
defined  earlier  and  employing  the  change  of  integration  variable  s  =  re±tn  to  calculate  the  function  4 !a(t )  = 
C~l{e~sl  °  },  t>  0,  where 

i  r°° 

(1.14)  #Q(t)  =  -  e~tre~r  '“«» ['(!-“)]  sin  (r1""  sin  [tt(1  -  a)])dr. 

ft  Jo 

Then,  using  (1.14),  a  change  of  variables  s'  =  s  in  (1.10)  and  subsequent  dropping  of  the  prime,  and 
the  translation  theorem  we  can  write  an  integral  expression  for  the  approximate  wavefront  response  (1.10) 
for  any  0  <  a  <  1  as  a  convolution  of  f(t)  with  the  approximate  impulse  response  function 

(1.15)  lapp{x,t)  =  C-1{e-As'~°}  =  A-^Va(A-T±;  (<-—)),  t  > 

Coo  Coo 

For  demonstration  purposes  we  employed  Mathematica  to  evaluate  the  integral  in  (1.14)  for  0  <  T  = 
t  -  -f-  <  2;  Figure  1.1  shows  the  wavefront  impulse  response  for  medium  parameters  and  depth  such  that 
A  =  1  and  various  values  of  a;  as  a  — >  1“  the  right  hand  side  of  (1.14)  approaches  a  delta  function  in  time 
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Fig.  1.3. 


Comparison  of  approximate  response  near  the  wavefront  to  a  numerical  solution  of  the  full  problem;  cx  T  =  0.036. 


and  the  numerical  evaluation  of  the  integral  fails.  We  note  that  in  the  case  of  unit  step-function  signaling 
data  (F(s)  =  1/s)  the  approximate  wavefront  response  exhibits  self-similarity,  i.e,  E(x,  t)  =  $«(£)>  where 
£  =  is  the  similarity  variable,  and  $Q(£)  —  Jq  l~a  t)dt,  where  we  have  used 

\pa(0)  =  0  and  the  integration  property  of  the  Laplace  transform.  Figure  1.2  shows  the  dependence  of 
lapp(x,t)  on  the  depth  x  (effected  by  varying  A)  for  a  =  0.7  and  representative  medium  parameters;  the 
rapidly  diminishing  vertical  scale  indicates  the  rapid  decay  of  the  wavefront  response.  We  also  note  that  the 
peak  of  the  response  behind  the  wavefront  T  =  0  is  delayed  by  an  amount  that  depends  on  both  a  and  x  (see 
Figures  1.1-1. 2).  We  have  verified  the  results  shown  in  Figures  1.1- 1.2  by  using  Mathematica  to  evaluate  the 
following  alternative  representation  of  # a(t )  [19], 

«■>« 

Figure  1.3  shows  a  validation  of  the  approximate  response  obtained  at  x  =  0.036cooT  by  convolving 
(1.15),  for  a  =  0.6,  with  a  rectangular  pulse  function  f(t)  of  duration  r  seconds.  The  numerical  result  (solid 
line)  is  a  solution  of  the  full  problem  consisting  of  (1.2)  and  (1.5)  using  a  novel  technique  we  developed  in 
[20]  for  computing  the  fractional  derivative  (1.6)  with  a  known,  a  priori  set,  error  that  is  uniform  in  time. 


The  t  «  ^  »  1  Region 

To  obtain  the  large-depth  impulse  response  we  evaluate  (1.8)  for  the  case  F(s)  =  1  using  the  saddle-point 
method.  We  first  rewrite  (1.8)  as 

i  rC+ioo 

(1.17)  E(x,t)  =  —  eXQ^ds, 

JQ-ioo 


where  6  =  >  1  is  the  space-time  parameter,  Q(s,0)  =  s[0  -  yj  J,  and  A  =  x/cQ 0  is  the  large 

parameter.  We  obtain  the  location  of  the  saddle  points  s  by  setting  a =  0  and  solving  for  s  =  s.  The 
following  equation  holds  for  the  saddle  points: 


(1.18) 


9  = 


sa  +  /3  _  a/3/2  a(3y/2 

sa+PY  sQ  +P  6‘a  +0Y 
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Fig.  1.4.  Motion  of  the  roots  of  (1.18)  for  l  <  0  <  ~^= 
and  a  =  0.7. 


Fig.  1.5.  Leading  order  asymptotic  response  at  x/coo 
100  as  a  function  of  a. 


By  rationalizing  (1.18)  we  obtain  a  fourth  order  polynomial  in  sa  whose  roots  we  find  numerically  for  a 
representative  set  of  parameters.  A  typical  root  portrait  is  shown  in  Figure  1.4;  the  arrows  indicate  the 
direction  of  motion  of  the  real  roots.  As  6  increases  from  1,  the  single  positive  real  root  sa  decreases  from  oo 
and  moves  towards  s  =  0  (the  arrow  labelled  with  the  positive  sign).  The  arrow  labelled  with  the  negative 
sign  indicates  the  direction  of  motion  of  the  negative  real  root  as  0  varies.  When  0  =  1/y/y  =  Cqo/cq  (as 
indicated  on  the  Figure),  the  real  root  is  sa  =  0.  For  0  >  Cqo/co  both  real  roots  of  (1.18)  are  negative.  The 
complex  roots  of  (1.18)  always  lie  in  the  second  and  third  quadrants.  Consequently,  on  the  principal  sheet 
defined  by  the  branch  cut  introduced  earlier  by  sa,  only  the  positive  real  root  of  (1.18)  survives  to  provide 
a  saddle  point.  Thus,  for  each  value  of  1  <  0  <  Cqo/co  there  corresponds  one  saddle  point  oo  >  s  >  0  on 
the  real  s— axis.  For  0  >  Cqq/cq  the  saddle  point  first  coalesces  with  the  branch  point  at  the  origin  and  then 
moves  into  the  branch  cut.  Consequently,  the  Bromwich  contour  is  no  longer  equivalent  to  a  steepest  descent 
contour.  In  the  region  1  <  6  <  Cqo/co  we  can  apply  the  saddle  point  method  using  A  =  x/cqo  as  our  large 
parameter.  A  small— s  expansion  of  Q(s,0), 


(1.19) 


allows  us  to  obtain  an  approximation  to  the  saddle  point,  i.e., 


(1.20) 


s  =  b«(4=-0)“> 


where  B  =  We  find  that  9  \s=s  =  &B  « (^  -  Q)1  i  >  0  hence  the  local  steepest  descent 

directions  at  s  are  arg(s  —  s)  =  ^ We  obtain  the  result 


(1.21) 


E(x,  t) 


pAQ(s,0) 


[1  + 


(l-a)(l  +  2a) 


y/27r\Qss(s,  0)  24aA Bl/*(-±*  -  6) 


i+i  J 


In  Figure  1.5  we  plot  the  leading  order  term  of  (1.21)  as  a  function  of  0  for  various  a,  and  note  that  its  peak 
does  not  occur  on  the  sub-characteristic  ray,  x  =  c$t  (0  =  ^),  as  it  does  in  the  case  of  the  Debye  medium 

(a  =  1)  [12]- [13].  Rather,  the  peak  of  the  response  arrives  earlier,  at  a  0  <  -^=,  and  its  space-time  location 
now  depends  on  a.  Also,  we  notice  that  the  leading-order  result  breaks  down,  i.e.,  E(x,t)  =  0,  for  9  =  -±= 
and  0  <  a  <  1;  the  second  term  in  (1.21),  the  correction,  diverges  at  6  =  -k=  since  the  derivation  of  (1.21) 
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Fig.  1.6.  Comparison  of  Equation  (1.21)  and  (1.1 7)  for  FlG.  1.7.  Comparison  of  Equation  (1.21)  and  (1.1 7)  for 

(a)  x/coo  =  50,  (b)  xj Coo  =  100.  The  arrow  indicates  the  (a)  x/coo  =  10,  (b)  x/coo  =  20,  (c)  x/coo  =  40.  .  The  arrow 
location  of  the  peak  of  the  response  for  a  =  1.  indicates  the  location  of  the  peak  of  the  response  for  a  =  1. 


Fig.  1.8.  Comparison  of  Equation  (1.21)  and  (1.1 7)  uhth  Fig.  1.9.  Large-depth  response  for  x/coo  =  18.432;  both 

a  numerical  simulation  [20]  for  x/coo  =  9.216.  The  arrow  in-  peak  location  and  value  depend  on  a.  The  arrow  indicates  the 

dicates  the  location  of  the  peak  of  the  response  for  a  =  1.  location  of  the  peak  of  the  response  for  a  =  1. 


does  not  take  into  account  the  coalescence  of  the  saddle  point  with  the  branch  point  at  s  =  0.  Despite  this 
failure,  the  result  is  useful;  Figures  1.6-1. 7  show  a  comparison  between  the  leading  order  sadle-point  method 
result,  and  an  evaluation  of  (1.17)  for  various  values  of  a  and  x/c0 0  using  Mathematica.  The  Mathematica 
results  have  been  verified  using  the  numerical  technique  developed  in  [20];  Figure  1.8  shows  such  a  verification 
where  the  numerical  result  overlaps  the  Mathematica  evaluation  of  (1.17),  as  indicated  by  the  boxed  region, 
and  both  results  agree  with  the  approximate  response  peak  location.  We  note  that  Figures  1.6- 1.8  indicate 
that  the  space- time  location  of  the  peak  response  also  depends  weakly  on  x/c0 G,  however  it  is  always  found 
in  the  region  6  <  ~^=. 

Although  we  have  not  been  able  to  approximate  the  response  shape  for  6  values  larger  than  those  corre¬ 
sponding  to  the  arrival  of  the  peak  response  at  a  given  depth,  we  have  investigated  this  region  numerically 
by  evaluating  (1.17)  with  Mathematica.  Figure  1.9  shows  our  results  for  0.5  <  a  <  0.9;  the  vertical  arrow 
indicates  the  arrival  of  the  peak  of  the  main  response  in  the  case  a  =  1  (Debye  model). 

1.1.4.  Summary.  We  have  determined  the  small-  and  large-depth  asymptotics  of  the  impulse  response 
of  the  Cole-Cole  dielectric  medium  model.  Our  results  were  validated  with  independent  numerical  solutions 
of  the  full  problem.  The  theoretical  wavefront  (short-depth)  and  large-depth  responses  of  the  Debye  and 
Cole- Cole  models  are  sufficiently  distinct  from  each  other  so  as  to  allow  propagated  pulse  measurements 
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in  TDR  setups  to  distinguish  between  the  two.  Significantly,  our  work  shows  that  the  measured  wavefront 
response  and  the  measured  arrival  of  the  main-response  peak  can  be  used  to  decide  which  dielectric  permit¬ 
tivity  model  is  most  appropriate  for  a  given  set  of  exerimentally-obtained  permittivity  data.  In  closing  we 
note  that  the  asymptotic  analysis  of  (1.17)  for  6  >  ^  remains  an  open  problem  which  must  be  solved  first 
in  order  to  obtain  a  uniform  asymptotic  approximation  of  the  impulse  response  for  1  <  6  <  oo. 


1.2.  Summary  of  Research  Results  to  be  Published. 

1.2.1.  The  failure  of  the  saddle-point  method.  The  asymptotic  analysis  of  (1.17)  was  shown  above 


from  below,  where  cQ 0  is  the  wavefront  speed 


to  be  accurate  for  1  <  6  <  ^  and  to  break  down  for  6 

—  Co 

and  Co  is  the  zero- frequency  phase  speed.  We  have  been  able  to  determine  the  asymptotic  form  of 

K+*oo 


(1.22) 


E(x,t)  =  -It  [  F(s)es|t"^V/i^lds,  f  >  — , 

2ni  7c -too  coo 


for  0  >  The  expression  we  obtained  is  as  follows: 


(1.23) 


where 


E(x,  t)  =  ^  — r(Q~^  Im{eiawF(  °  +  1 


a+1 


x(0 


-J£' 

7  • 


,«)}, 


A  i /£(!_!), 

2y  1  P  r 


and 


F(r,  a)  - 


rAeio 


We  are  presently  attempting  to  verify  (1.23)  by  comparing  with  numerical  simulations  obtained  with  the 
method  described  below  and  to  fill  the  gap  in  the  analysis  for  9  =  ^ . 

1.2.2.  Numerical  methods  for  fractional  differential  equations.  To  incorporate,  e.g.,  the  Cole- 
Cole  model  in  the  FD-TD  scheme  we  must  be  able  to  numerically  solve  the  following  fractional  differential 
equation  initial  value  problem  for  P 

P 


+  P  =  (e,  -  eoo)E;  t  >  0,  P(0)  =  0 


where 


da p  _  i  rl  p'(i) 
dta  r(l  -  a)  J0  (t  -  £)a 

is  the  Caputo  fractional  derivative  of  order  a  €  (0, 1).  We  use 


di 


ro 


«)=  f° 

Jo 


e~zza~1 


dz  and  r(a)r(l  —  a)  =  — 


Sin  7TQ 


and  a  simple  change  of  variables,  z  =  £2(t  -  £),  in  the  Caputo  derivative  to  show: 


a=0.7 


(a)  (b) 

Fig.  1.10.  (a)  Demonstration  of  accuracy  of  scheme  to  compute  fractional  derivatives;  (b)  demonstration  of  accuracy 
obtained  with  the  scheme  near  t  =  0. 

where  t/>(£,  t)  satisfies  a  non-homogeneous  ordinary  differential  equation  initial  value  problem 

^+£V  =  S2“-1^;  <>o,  ip(£,o)  =  0;  o<^<^. 

For  accurate  results  we  decompose  J0°°  d£  =  f*  d£  +  and  use  the  small— £  behavior  i/>(£,  t)  ~  y2a~lP(t) 
and  the  large-£  behavior  ^>(£,  t)  ~  y2a~3P'  (t),  determined  from  the  above  differential  equation,  along  with 
Gaussian  and  Laguerre  integration  respectivelly.  The  discretization  in  £  is  then  followed  by  a  discretization 
in  t  using  the  trapezoidal  rule;  the  result  is  an  0(At2)  accurate  system  of  Ng  4-  Ni  4- 1  equations  that  must 
be  solved  along  with  the  FD-TD  scheme.  Typically,  Ng  =  Ni  =  10  -  15.  The  validation  of  the  approach  is 
shown  in  Figure  1.10  where  the  solution  of  the  fractional  ordinary  differential  equation 

^  +  P-0,  P<0+)  =  1, 

is  compared  to  the  exact  solution 

P (t)  =  Ea(—ta),  t  >  0. 

This  work  is  in  review  [20]. 

1.2.3.  The  Cole-Davidson  and  Havriliak-Negami  dielectric  models.  Previously,  we  found  that 
the  Debye,  €d(u)  =  ,  wavefront  impulse  response  supports  discontinuities  on  x  =  ,  where 

Coo  =  l/\A(°o )mo  is  the  wavefront  speed,  and  that  it  decays  exponentially  with  depth  hence  it  is  important 
only  up  to  a  distance  of  O(coot)  past  the  air /dielectric  interface  (the  so-called  ’’time-domain  skin-depth”). 
Also,  we  found  that  past  that  thin  layer  the  impulse  response  satisfies  an  advection-diffusion  equation, 
that  it  travels  with  the  zero- frequency  phase  speed  along  the  sub-characteristic  x  =  =  cot,  where 

es  =  6(0)  is  the  DC  permittivity,  and  that  it  diffuses  around  that  sub-characteristic.  The  Cole-Cole  model’s 
impulse  response  is  infinitely  smooth  at  the  wavefront  (small-depth).  This  is  in  contrast  to  the  case  of 
the  Debye  impulse  response  which  is  discontinuous  at  the  wavefront.  Also,  the  location  of  the  peak  of  the 
main  (large-depth)  response  occurs  at  an  earlier  space-time  location  than  that  found  for  the  Debye  dielectric 
model  main  response.  Hence,  the  short-depth  and  large-depth  responses  of  the  Debye  and  Cole-Cole  models 
are  sufficiently  distinct  from  each  other  so  as  to  allow  propagated  pulse  measurements  in  TDR  setups  to 
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distinguish  between  the  two  when  used  to  numerically  model  such  experiments  by  fitting  the  two  models  to 
the  same,  experimentally  determined,  dielectric  permittivity  data. 

For  this  reason  we  decided  to  extend  the  work  done  on  the  Cole-Cole  dielectric  model  in  order  to 
investigate  the  remaining  dielectric  models  of  anomalous  dispersion,  namely  the  Cole-Davidson  and  Havriliak- 
Negami  models.  The  Maxwell  system  in  these  dielectrics  takes  the  following  form 

Df  =  V  x  H,  /ioHt  =  -V  x  E 
and  is  closed  with  the  time-domain  contitutive  law 

D(£)  =  CooE(£)  +  [  ^  )E(*  )dt 

Jo 


where  (a  G  (0, 1),  (3  G  (0, 1)) 


X(t)  =  XCD(t)  =  t  >  0 

x(t)  =  xHN(t)  =  1  >  0 


or 


x(t)  =  xD(t)  =  c-1C- <>0- 

1  4-  .ST 

Therefore  we  undertook  to  examine  wave  propagation  (and  develop  numerical  schemes  to  solve  for  the 
propagation  of  EM  signals)  in  the  following  models: 

•  ecc(^)  =  too  4-  ,  0  <  a  <  1,  was  proposed  by  Cole  &  Cole  (J.  Chem.  Phys.,  v.  9,  1941)  as 

an  alternative  to  the  Debye  model  (a  =  1). 

•  ccd(^)  =  foo  4-  (ii7J?)34  0</?<l,  was  proposed  as  an  alternative  to  the  above  model  by  Davidson 
&  Cole  (J.  Chem.  Phys.,  v.  19,  1951)  in  order  to  have  a  model  that  resembles  the  Debye  medium 
in  the  low-frequency  limit  while  providing  for  a  distribution  of  relaxation  times. 

•  chn{v)  =e oo  4-  »  0  <  a  <  1,  0  <  /?  <  1,  was  proposed  by  Havriliak  &  Negami  (J.  Polym. 

Sci.,  v.  14,  1966)  to  combine  the  features  in  the  above  two  models. 

Alternatively,  the  constitutive  law  above  can  be  implemented  via  writing 

D(t)  =  €ooE(t)  -f  P (t) 

where  P (t)  (P(0)  =  0)  satisfies  a  fractional  differential  equation  that  involves 

da P  =  1  P \t) 

dta  T(1  -  a)  J0  (t  -  0Q 

which  is  the  Caputo  fractional  derivative  of  order  a  G  (0, 1).  For  numerical  purposes  one  has  to  numerically 
solve  fractional  differential  equations. 
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First,  some  results  for  the  Cole-Cole  Susceptibility  obtained  in  our  previous  work.  A  representation  of 
the  Cole-Cole  susceptibility  is 


xccw  =  A]T(-ir+i 

n=  1 


T(na) 


t  >  0 


which,  by  referring  to  textbooks,  can  be  written  as 

Xcc(t)  =  -TAftEa(-(t/T)a),  t  >  0 

where  Ea(z)  =  r(i+ noTj  *s  the  Mittag-Leffler  function  of  order  a.  For  example,  using  a  table  of 

transforms,  when  a  =  1/2 

XCC(t)  =  A[  J—  -  et^Ter/c( \ftfr )] ,  t  >  0 
yirf/r 

Given  a,  t,  and  r  one  finds  that  O(^)  terms  must  be  summed  in  order  to  evaluate  Xa(t)  and  perform  the 
convolution  numerically  in  a  wave  propagation  code.  In  addition,  for  0  <  a  <  1,  we  have 

xcc(t)~o(^),  t-0+ 


XCC(t)  ~  O(^L),  t  —  00 

Thus,  xCC(t)  €  Lp  for  1  <  p  <  thus  only  for  a  >  1/2  does  x<*(t)  belong  to  the  intersection  of  L1  and 
L2,  i.e.,  the  Fourier  transform  (and  the  Kramers-Kronig  relations)  require  more  care  for  a  <  1/2  since  then 
XCC(t)  is  only  in  L1.  In  biological  and  geophysical  applications  0.6  <  a  <  0.9  . 

Next,  we  give  some  results  for  the  Cole-Davidson  Susceptibility  model.  For  the  Cole-Davidson  dielectric, 
the  time-domain  susceptibility  is  (0  — >  1“  is  the  Debye  model) 


<CD(t)  =  c-\- 


)  = 


t**- 


k(l  +  st)P 

which  can  be  used  to  determine  the  polarization  ODE 


m 


e~l/T ,  t  >  0 


^  +  ^p  =  e 


■t/T^[et/Tp]  =  t  >  0,  P(0)  =  0 


for  use  in  the  constitutive  law 


D(0=cooE(t)  +  P(t) 

where  ^  is  again  the  Caputo  fractional  derivative  of  order  0  £  (0, 1).  The  Cole-Davidson  fractional 
differential  equation  was  obtained  through  use  of  the  following  Laplace  transform  identity 

A  A@ 

£{(!  +  = 


C{^(e^f(t))}(s  +  7)  =  (s  +  7 +  7)  -  /( 0+)(s  +  7)*-1  = 
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{s  +  ifFis), 

since  our  initial  condition  on  the  polarization  requires  /(0+)  =  0,  and  the  following  relation 

es  -  Coo  E(s) 


PW  = 


(1  4-  st)P 


for  the  polarization  in  the  frequency  domain.  Numerically,  the  Cole- Davidson  model  would  employ  the 
following 


where  ipCD(£,t)  satisfies  a  non-homogeneous  ordinary  differential  equation  initial  value  problem 

diPCD 


dt 


+  fV 


2,/, CD  _  (10-\ ipul/r. 


-  +  -P)e*'r;  t>  0, 
at  t 


subject  to  the  initial  condition 


*0  (^,  0)  =  0;  0  <  £  <  oo 


in  order  to  solve 

^[e‘/TP]  =  ££-^Eet/T;  t  >  0,  P(0)  =  0 
In  order  to  tackle  the  Havriliak-Negami  model  we  need  the  following: 


i  (-^H) 

rfc  r(i  +  fc)r(— /?) 


so  that  the  application  of  the  inverse  Laplace  transform  to  both  sides  gives 


,±  ,  Iw»  =  V  1  (-i)fcr(fc-/9)  d  fl_fc 

''dt  r;  ^'|rfcr(l  +  A:)r(-^)''dt; 


Then  we  can  solve  the  fractional  differential  equation  (H(t)  is  the  Heavyside  step  function) 

<S  +  ;  «')-> 

using  the  Laplace  transform  and  subsequently  obtain  the  Green’s  function,  ,  for  solving  it  with  a  general 
right-hand  side,  i.e., 


A  (~l)k+l(P  +  k)T(P  -f  k)  t  ,0+k-i 

T  F(/?)r(A;  4-  l)r(l  +  ft  +  ky  r 


=  X 


CD 


<‘>  =  7£ 


fc=0 


(_i)fc+i  *  j3+fc_1 

r(/?)r(fc  +  l)  V 


which  is  just  an  alternative  representation  of  ^Y(p)e  •  Clearly, 
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00 


XCD{t)  ~  0{e~t/T),  t- 


So  the  CD  model  resembles  the  CC  at  short  times  and  the  Debye  at  long  times.  In  the  above  derivation  we 
used  £;H(t)  = 

Finally,  we  give  results  obtained  for  the  Havriliak-Negami  susceptibility  function.  To  find  the  Green’s 
function  for  the  fractional  differential  equation  which  will  be  coupled  to  Maxwell’s  equations  we  use 


oo 


k=0 


_J_  (~l)*T(fc  +  /?)  a(0+k) 
rka  T(l+k)T{0) 


and  the  property  ^ H(t )  =  J  ^ ,  to  first  solve 


p  W  =  aE 


k= 0 


(— i)fc+1r(/?  +  fc)  t  (0+k) 
r(/?)r(it  +  i)r(i  +  Q(/3  +  fc)V; 


Subsequent  time  differentiation  gives 

vHN(t)  =  —  V'  (— 1)fc+1r'(/3  +  k)  (l^a(0+k)-i 

T  ~r(flr(k  +  i)r(a(/9  +  *)V 

The  asymptotic  behavior  now  is 

xHN(t)~0( t-0+ 


X™(t)~0(^).  t-oo 

So  the  HN  model  resembles  the  CC  at  short  times  (but  with  slower  decay)  and  at  long  times  (identical 
decay).  For  the  non-dimensionalized  signaling  problem  in  a  Havriliak-Negami  dielectric  half  space  (x  >  0) 
the  electric  field  is  given  by  the  following  Bromwich  integral  (E  =  0,  t  <  x) 


where  -7r  <  arg(s)  <  7r,  fie/  >  0,  7  =  (^)a/3(^  -  1),  6  =  and  F(s)  is  the  Laplace  transform  of 
the  time-domain  signaling  data  imposed  at  x  =  0.  For  the  branch  chosen,  the  SQRT  does  not  introduce 
additional  singularities  (recall,  0<a<l&0</3<l).  We  first  want  to  know  the  t  — ►  x+  behavior  of  the 
response  with  F(s)  =  1,  i.e.,  the  wavefront  behavior  of  the  HN  Green’s  function.  The  large— s  expansion 


t-x  +  x(  1  -  /+  (g0  //)  =  (t-x)~  x^s  a0  +  0(xs  Q(1+/J)) 


gives 


-As1-00 
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ea(t~x)ds, 


where  A  =  >  0.  Since  1  —  ol0  >  0,  7  >  0,  we  find 

dn 

lim  - — E{x,t )  =0,  n  >  0, 

*— x+  dtn 

since 

lim  snF(s)e~Al>l  al3  =  0,  n  >  0. 

s  — >oo 

Thus,  the  response  of  the  Havriliak-Negami  medium  is  infinitely  smooth  at  the  wavefront  even  in  the  case  of 
F(s)  =  1  ( f(t )  =  S(t)  signaling  data).  This  is  similar  to  the  Cole-Cole  dielectric.  In  contrast,  for  the  Debye 
model  (a  =  1,  0  =  1)  we  obtain 

E(x,t )  ~  e~  2c^r(<7^~1^  f(t - — ), 

Cqo 

i.e. ,  the  Debye  wavefront  response  inherits  the  continuity  properties  of  f(t)  at  t  =  0  and  decays  exponentially 
fast  thus  exhibiting  a  ’’time-domain  skin-depth”  (see  TMR  &  PGP,  JOSA-A  (1996),  and  PGP,  WAVE 
MOTION  (1995)).  Integrating  C~l{e~sX  a0}  around  the  branch  cut  we  obtain  the  wavefront  response  for 
any  a,  0  in  terms  of  the  function 

i  r°° 

Va0{T)  =  -  I  e~Tre~r  a0  cos  sin  (rl~a(3  sin  [7r(l  —  a0)])dr, 

n  Jo 

T  >  0,  where  T  =  t  —  x.  NOTE:  lima<  ^fa,p{T)  =  e~18(T).  Then,  the  translation  theorem  allows  us 
to  write  an  integral  expression  for  the  approximate  wavefront  response  for  any  0  <  a  <  1,  0  <  /?  <  1  as  a 
convolution  of  f(t)  with  the  approximate  impulse  response  function 

lapp{x,t)  =  £-1{e"iUl'-'}  =  A-r^<nnJ3{A-r^(t  -  x)),  t  >  x. 

A  publication  that  compares  and  contrasts  these  various  dielectric  models  from  the  point  of  view  of  their 
respective  EM  pulse  responses  is  in  preparation. 
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